|--------------------------------------------------|
|==================================================|
'select()' returned 1:many mapping between keys and columns
'select()' returned 1:many mapping between keys and columns
'select()' returned many:many mapping between keys and columns
#QC
read.depth.plots(rnaseq$meta)
Slightly higher mouse fraction than normal, but overall, this is damn near perfect.
pca <- pca.run(rnaseq$human)
pca <- cbind(rnaseq$meta,pca)
plts <- list(
ggpairs(pca, columns=c('PC1', 'PC2', 'PC3', "human.library.size"),mapping=aes(col=Treatment.arm), upper=list(continuous='cor')),
ggpairs(pca, columns=c('PC1', 'PC2', 'PC3', "human.library.size"),mapping=aes(col=Group), upper=list(continuous='cor'))
)
for(p in plts) print(p, progress = F)
The largest differences are between the experimental groups (short treatment, long treatment, & resistant samples), with the application of the drug having a smaller, if any, effect. What do the PCAs look like if we look at each group inidividually? Also what are the outliers for group 1?
plts <- lapply(levels(rnaseq$meta$Group),function(g){
rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Group == g])
pca <- pca.run(rnaseq.subset$human)
pca <- cbind(rnaseq.subset$meta,pca)
ggpairs(pca, columns=c('PC1', 'PC2', 'PC3', "human.library.size"),mapping=aes(col=Treatment.arm), upper=list(continuous='cor'),title = paste("Group:",g))
})
for(p in plts) print(p, progress = F)
There appears to be an outlier in group 1, which we will identify below. Otherwise, nothing of note, which indicates that there may not be large transcriptomic changes.
rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Group == 1])
pca <- pca.run(rnaseq.subset$human)
pca <- cbind(rnaseq.subset$meta,pca)
pca
Sample.name Pool Barcode Sequence Flowcell Lane fname mouse.library.size human.library.size Group Origin Type.of.material Treatment.arm
1: HCI011-T1-P01-18.241B-T1 SLX-16182 UDI0065 ACACTAAG-ATCCATAT H23KHBBXY s_4 SLX-16182.UDI0065.H23KHBBXY.s_4 726453 10740675 1 HCI011-x3 (spare of tumour of HCI011-x2 implanted) 3 doses vehicle
2: HCI011-T1-P01-18.9413L2-T1 SLX-16182 UDI0066 GTGTCGGA-GCTTGCGC H23KHBBXY s_4 SLX-16182.UDI0066.H23KHBBXY.s_4 598305 9716980 1 HCI011-x2 3 doses vehicle
3: HCI011-T1-P01-18.9413R1-T1 SLX-16182 UDI0067 TTCCTGTT-AGTATCTT H23KHBBXY s_4 SLX-16182.UDI0067.H23KHBBXY.s_4 266097 12349091 1 HCI011-x2 3 doses vehicle
4: HCI011-T1-P01-18.9416L2-T1 SLX-16182 UDI0068 CCTTCACC-GACGCTCC H23KHBBXY s_4 SLX-16182.UDI0068.H23KHBBXY.s_4 1199394 10065053 1 HCI011-x2 3 doses vehicle
5: HCI011-T1-P01-18.239L1-T1 SLX-16182 UDI0069 GCCACAGG-CATGCCAT H23KHBBXY s_4 SLX-16182.UDI0069.H23KHBBXY.s_4 1056432 9397878 1 HCI011-x3 (spare of tumour of HCI011-x2 implanted) 3 doses GDC0032
6: HCI011-T1-P01-18.242B-T1 SLX-16182 UDI0070 ATTGTGAA-TGCATTGC H23KHBBXY s_4 SLX-16182.UDI0070.H23KHBBXY.s_4 1184487 9929416 1 HCI011-x3 (spare of tumour of HCI011-x2 implanted) 3 doses GDC0032
7: HCI011-T1-P01-18.242L2-T1 SLX-16182 UDI0071 ACTCGTGT-ATTGGAAC H23KHBBXY s_4 SLX-16182.UDI0071.H23KHBBXY.s_4 1574978 10833734 1 HCI011-x3 (spare of tumour of HCI011-x2 implanted) 3 doses GDC0032
8: HCI011-T1-P01-18.242NM-T1 SLX-16182 UDI0072 GTCTACAC-GCCAAGGT H23KHBBXY s_4 SLX-16182.UDI0072.H23KHBBXY.s_4 833155 10428969 1 HCI011-x3 (spare of tumour of HCI011-x2 implanted) 3 doses GDC0032
Days.treatment Sensitivity Volume.of.Tumor Cause.of.death PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
1: 4 - 498 end of exp -4.865990 11.206267 3.3524440 -3.6778368 6.7825653 -2.87223053 5.2110471 2.220446e-16
2: 4 - 695 end of exp -3.900719 9.939373 14.8538190 -1.3780865 -4.4327113 2.07627946 -1.5520391 6.383782e-16
3: 4 - 540 end of exp 59.960192 -6.246141 -0.9891759 0.2395942 -0.2123211 0.05816152 0.1011068 -2.579881e-14
4: 4 - 1066 end of exp -8.522011 8.161411 -3.3656177 7.4147375 -5.0739557 -5.75831411 0.8230278 1.600109e-14
5: 4 Sensitive 570 end of exp -17.093052 -37.367157 2.9208371 2.1070527 1.4118816 0.14123286 0.5240586 4.898859e-15
6: 4 Sensitive 464 end of exp -10.586991 -2.650520 -8.1898604 -9.9072927 -4.8273185 -0.14422147 -0.3452060 8.978929e-15
7: 4 Sensitive 648 end of exp -8.577357 9.543896 -6.0345708 4.6027311 0.4754336 7.33337173 2.3050118 4.662937e-15
8: 4 Sensitive 498 end of exp -6.414070 7.412872 -2.5478753 0.5991005 5.8764260 -0.83427945 -7.0670070 3.802514e-15
y <- pca.prep(rnaseq.subset$human,T,1000)
y <- prcomp(t(y),scale=T)
summary(y)
Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8
Standard deviation 24.5698 16.3854 7.19815 5.2726 4.65214 3.78413 3.50308 1.02e-14
Proportion of Variance 0.6037 0.2685 0.05181 0.0278 0.02164 0.01432 0.01227 0.00e+00
Cumulative Proportion 0.6037 0.8721 0.92397 0.9518 0.97341 0.98773 1.00000 1.00e+00
Outlier is “HCI011-T1-P01-18.9413R1-T1” w.r.t to PC1, whose only notable difference is a very low fraction of mouse (~2%).
plts <- lapply(unique(rnaseq$meta$Treatment.arm),function(ta){
rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Treatment.arm == ta])
pca <- pca.run(rnaseq.subset$human)
pca <- cbind(rnaseq.subset$meta,pca)
ggpairs(pca, columns=c('PC1', 'PC2', 'PC3', "human.library.size"),mapping=aes(col=Group), upper=list(continuous='cor'),title = paste("Treatment arm:",ta))
})
for(p in plts) print(p, progress = F)
rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Group == 1])
#rnaseq.subset <- subset.PDX(rnaseq.subset, rnaseq.subset$meta[,Sample.name!="HCI011-T1-P01-18.9413R1-T1"])
treatment <- rnaseq.subset$meta[,factor(Treatment.arm,levels=c("vehicle","GDC0032"))]
d <- model.matrix(~ 0 + treatment)
d
treatmentvehicle treatmentGDC0032
1 1 0
2 1 0
3 1 0
4 1 0
5 0 1
6 0 1
7 0 1
8 0 1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment
[1] "contr.treatment"
cont.matrix <- matrix(c(-1,1), dim(d)[2], 1)
colnames(cont.matrix) <- c("GDC.minus.control")
rownames(cont.matrix) <- colnames(d)
cont.matrix
GDC.minus.control
treatmentvehicle -1
treatmentGDC0032 1
#Perform DE
DE.data <- DE.run(rnaseq.subset$human,d,cont.matrix = cont.matrix)
#Top DE genes table
DE <- topTable.annotated.DE(DE.data$efit,p.val = 1)
'select()' returned 1:many mapping between keys and columns
if(nrow(DE[adj.P.Val<0.05])!=0){
DE <- DE[adj.P.Val<0.05]
table(DE[,direction])
DE[direction=="up.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
DE[direction=="down.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
#Pheatmap w/ batch correction
i <- match(DE$ENSEMBL,rownames(DE.data$v$E))
y <- removeBatchEffect(DE.data$v$E[i,],design = model.matrix(~drug), batch= pdx)
labels.col <- rep("",nrow(rnaseq$meta))
labels.row <- DE$SYMBOL
ann.col <- as.data.frame(rnaseq$meta[,.(pdx,drug)])
rownames(ann.col) <- rnaseq$meta[,Sample.name]
pheatmap(y,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
} else {
head(DE,10)
}
ENSEMBL logFC AveExpr t P.Value adj.P.Val B direction ENTREZID SYMBOL GENENAME
1: ENSG00000204388 1.0760352 9.253658 4.832736 0.0006147231 0.3272351 -1.519333 up.reg 3304 HSPA1B heat shock protein family A (Hsp70) member 1B
2: ENSG00000236675 -0.8133502 5.721730 -4.733665 0.0007165805 0.3272351 -1.574151 down.reg <NA> <NA> <NA>
3: ENSG00000183011 -0.8991874 6.269441 -4.568031 0.0009288204 0.3272351 -1.620282 down.reg 84316 NAA38 N(alpha)-acetyltransferase 38, NatC auxiliary subunit
4: ENSG00000145494 -0.7980749 6.345503 -4.512080 0.0010147699 0.3272351 -1.652281 down.reg 4726 NDUFS6 NADH:ubiquinone oxidoreductase subunit S6
5: ENSG00000237991 -1.0496058 6.159797 -4.452303 0.0011159358 0.3272351 -1.719579 down.reg <NA> <NA> <NA>
6: ENSG00000151208 1.0660589 5.918073 4.490701 0.0010497920 0.3272351 -1.728342 up.reg 9231 DLG5 discs large MAGUK scaffold protein 5
7: ENSG00000077809 0.8711585 6.583563 4.340502 0.0013347542 0.3272351 -1.771400 up.reg <NA> <NA> <NA>
8: ENSG00000105669 -0.7967949 7.069990 -4.264799 0.0015082447 0.3272351 -1.807395 down.reg 11316 COPE coatomer protein complex subunit epsilon
9: ENSG00000125356 -0.6583621 6.864270 -4.197563 0.0016822276 0.3272351 -1.863827 down.reg 4694 NDUFA1 NADH:ubiquinone oxidoreductase subunit A1
10: ENSG00000180389 -0.7912983 6.400057 -4.216161 0.0016320902 0.3272351 -1.873863 down.reg 432369 ATP5F1EP2 ATP synthase F1 subunit epsilon pseudogene 2
#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- merge(y,rnaseq.subset$meta)
ggplot(y) + aes(x=SYMBOL,y=Expression,color=Treatment.arm) + geom_boxplot() + stat_compare_means(aes(label=paste0("p.adj = ",..p.adj..)),method="t.test") + theme_tufte() + theme(legend.position = "top")
## Heatmap
labels.col <- rnaseq.subset$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq.subset$meta[,.(Treatment.arm)])
rownames(ann.col) <- rnaseq.subset$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
#GSEA
## C6
gsea <- GSEA.run(DE.data$v,d,cont.matrix[,1],Hs.c6,0.05)
'select()' returned 1:many mapping between keys and columns
gsea[Direction=="Up",.(rn,NGenes,Direction,FDR)]
rn NGenes Direction FDR
1: TGFB_UP.V1_UP 133 Up 5.785652e-05
2: STK33_NOMO_UP 210 Up 2.290397e-04
3: STK33_UP 201 Up 3.175661e-04
4: LTE2_UP.V1_DN 174 Up 3.972324e-04
5: STK33_SKM_UP 192 Up 9.710797e-04
6: BCAT_BILD_ET_AL_DN 38 Up 9.766632e-04
7: EIF4E_DN 90 Up 2.000560e-03
8: BMI1_DN.V1_UP 101 Up 2.347305e-03
9: BMI1_DN_MEL18_DN.V1_UP 95 Up 2.444336e-03
10: TBK1.DF_DN 266 Up 5.350087e-03
11: CAMP_UP.V1_DN 155 Up 5.671311e-03
12: CSR_EARLY_UP.V1_DN 121 Up 6.646370e-03
13: PIGF_UP.V1_UP 152 Up 8.964812e-03
14: EGFR_UP.V1_UP 168 Up 8.964812e-03
15: MEK_UP.V1_DN 172 Up 1.038029e-02
16: RAF_UP.V1_DN 177 Up 1.038029e-02
17: MEL18_DN.V1_UP 96 Up 1.299166e-02
18: RPS14_DN.V1_UP 97 Up 1.354666e-02
19: PRC2_EZH2_UP.V1_DN 114 Up 1.940560e-02
20: PDGF_UP.V1_DN 70 Up 1.940560e-02
21: PDGF_ERK_DN.V1_DN 130 Up 2.272031e-02
22: RAF_UP.V1_UP 167 Up 2.272031e-02
23: RB_P130_DN.V1_DN 121 Up 2.272031e-02
24: KRAS.DF.V1_DN 129 Up 2.272031e-02
25: HINATA_NFKB_IMMU_INF 7 Up 2.297676e-02
26: MTOR_UP.N4.V1_DN 146 Up 2.652461e-02
27: PDGF_UP.V1_UP 127 Up 3.046860e-02
28: E2F1_UP.V1_DN 148 Up 3.196333e-02
29: HOXA9_DN.V1_UP 148 Up 3.495692e-02
30: PRC1_BMI_UP.V1_DN 86 Up 3.584935e-02
31: PRC2_EED_UP.V1_UP 108 Up 3.584935e-02
32: ALK_DN.V1_DN 64 Up 4.019820e-02
33: RELA_DN.V1_UP 92 Up 4.672090e-02
34: MEK_UP.V1_UP 170 Up 4.672090e-02
35: IL21_UP.V1_DN 74 Up 4.672090e-02
36: ERB2_UP.V1_DN 185 Up 4.672090e-02
37: MYC_UP.V1_DN 116 Up 4.672090e-02
38: SINGH_KRAS_DEPENDENCY_SIGNATURE_ 18 Up 4.672090e-02
rn NGenes Direction FDR
gsea[Direction=="Down",.(rn,NGenes,Direction,FDR)]
rn NGenes Direction FDR
1: EIF4E_UP 77 Down 0.0002290397
2: CSR_LATE_UP.V1_UP 142 Down 0.0002463503
3: RB_P107_DN.V1_UP 116 Down 0.0002496484
4: RB_P130_DN.V1_UP 100 Down 0.0194056047
5: MYC_UP.V1_UP 142 Down 0.0227203071
6: RPS14_DN.V1_DN 155 Down 0.0236125612
rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Group == 2])
treatment <- rnaseq.subset$meta[,factor(Treatment.arm,levels=c("vehicle","GDC0032"))]
d <- model.matrix(~ 0 + treatment)
d
treatmentvehicle treatmentGDC0032
1 1 0
2 1 0
3 1 0
4 0 1
5 0 1
6 0 1
7 0 1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment
[1] "contr.treatment"
cont.matrix <- matrix(c(-1,1), dim(d)[2], 1)
colnames(cont.matrix) <- c("GDC.minus.control")
rownames(cont.matrix) <- colnames(d)
cont.matrix
GDC.minus.control
treatmentvehicle -1
treatmentGDC0032 1
#Perform DE
DE.data <- DE.run(rnaseq.subset$human,d,3,cont.matrix = cont.matrix)
#Top DE genes table
DE <- topTable.annotated.DE(DE.data$efit,p.val = 1)
'select()' returned 1:many mapping between keys and columns
if(nrow(DE[adj.P.Val<0.05])!=0){
DE <- DE[adj.P.Val<0.05]
table(DE[,direction])
DE[direction=="up.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
DE[direction=="down.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
#Pheatmap w/ batch correction
i <- match(DE$ENSEMBL,rownames(DE.data$v$E))
y <- removeBatchEffect(DE.data$v$E[i,],design = model.matrix(~drug), batch= pdx)
labels.col <- rep("",nrow(rnaseq$meta))
labels.row <- DE$SYMBOL
ann.col <- as.data.frame(rnaseq$meta[,.(pdx,drug)])
rownames(ann.col) <- rnaseq$meta[,Sample.name]
pheatmap(y,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
} else {
print("No DE genes found")
head(DE,10)
}
[1] "No DE genes found"
ENSEMBL logFC AveExpr t P.Value adj.P.Val B direction ENTREZID SYMBOL GENENAME
1: ENSG00000174021 1.0508825 6.847082 6.713034 7.674007e-05 0.6224003 -3.671849 up.reg 2787 GNG5 G protein subunit gamma 5
2: ENSG00000164128 -1.1130122 5.649225 -5.955737 1.919213e-04 0.9268516 -3.790498 down.reg 4886 NPY1R neuropeptide Y receptor Y1
3: ENSG00000204569 0.8046126 6.989129 4.701910 1.037803e-03 0.9268516 -3.835360 up.reg 5514 PPP1R10 protein phosphatase 1 regulatory subunit 10
4: ENSG00000164442 -0.9972152 5.227045 -5.581024 3.105660e-04 0.9268516 -3.852576 down.reg 10370 CITED2 Cbp/p300 interacting transactivator with Glu/Asp rich carboxy-terminal domain 2
5: ENSG00000186918 -0.6304149 5.900330 -4.502097 1.386601e-03 0.9268516 -3.907247 down.reg 55893 ZNF395 zinc finger protein 395
6: ENSG00000104267 -0.8459846 6.814805 -4.081274 2.600769e-03 0.9268516 -3.919892 down.reg 760 CA2 carbonic anhydrase 2
7: ENSG00000144802 -1.1270393 5.645451 -4.275787 1.938616e-03 0.9268516 -3.952431 down.reg 64332 NFKBIZ NFKB inhibitor zeta
8: ENSG00000074800 -0.4960966 8.821536 -3.652362 5.062147e-03 0.9268516 -3.974104 down.reg 2023 ENO1 enolase 1
9: ENSG00000135127 0.8043419 6.155308 3.951610 3.172693e-03 0.9268516 -3.976310 up.reg 92558 BICDL1 BICD family like cargo adaptor 1
10: ENSG00000104321 -1.8460421 7.281383 -3.667565 4.942057e-03 0.9268516 -3.983110 down.reg 8989 TRPA1 transient receptor potential cation channel subfamily A member 1
#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- merge(y,rnaseq.subset$meta)
ggplot(y) + aes(x=SYMBOL,y=Expression,color=Treatment.arm) + geom_boxplot() + stat_compare_means(aes(label=paste0("p.adj = ",..p.adj..)),method="t.test") + theme_tufte() + theme(legend.position = "top")
## Heatmap
labels.col <- rnaseq.subset$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq.subset$meta[,.(Treatment.arm)])
rownames(ann.col) <- rnaseq.subset$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
#GSEA - c6
gsea <- GSEA.run(DE.data$v,d,cont.matrix[,1],Hs.c6,0.05)
'select()' returned 1:many mapping between keys and columns
gsea[Direction=="Up",.(rn,NGenes,Direction,FDR)]
rn NGenes Direction FDR
1: ESC_V6.5_UP_EARLY.V1_DN 116 Up 8.734097e-06
2: BMI1_DN.V1_UP 101 Up 1.146764e-03
3: BMI1_DN_MEL18_DN.V1_UP 94 Up 1.299194e-03
4: TGFB_UP.V1_UP 133 Up 3.186118e-03
5: KRAS.600_UP.V1_DN 134 Up 3.186118e-03
6: PTEN_DN.V2_UP 78 Up 4.432734e-03
7: KRAS.AMP.LUNG_UP.V1_DN 74 Up 5.805192e-03
8: KRAS.300_UP.V1_DN 72 Up 8.586000e-03
9: KRAS.LUNG_UP.V1_DN 68 Up 1.156810e-02
10: MEL18_DN.V1_UP 95 Up 1.211232e-02
11: SNF5_DN.V1_DN 92 Up 3.277517e-02
12: KRAS.KIDNEY_UP.V1_DN 60 Up 3.277517e-02
13: CAHOY_ASTROGLIAL 59 Up 3.277517e-02
14: RELA_DN.V1_DN 58 Up 3.277517e-02
15: MEK_UP.V1_UP 171 Up 3.277517e-02
16: CYCLIN_D1_UP.V1_UP 127 Up 3.874861e-02
17: CAHOY_OLIGODENDROCUTIC 69 Up 4.393997e-02
18: P53_DN.V1_UP 147 Up 4.393997e-02
19: SINGH_KRAS_DEPENDENCY_SIGNATURE_ 18 Up 4.393997e-02
20: RAPA_EARLY_UP.V1_UP 115 Up 4.393997e-02
21: MYC_UP.V1_DN 115 Up 4.589087e-02
22: VEGF_A_UP.V1_UP 129 Up 4.589087e-02
23: AKT_UP.V1_DN 132 Up 4.589087e-02
rn NGenes Direction FDR
gsea[Direction=="Down",.(rn,NGenes,Direction,FDR)]
rn NGenes Direction FDR
1: HOXA9_DN.V1_DN 160 Down 0.03277517
2: TBK1.DF_DN 268 Down 0.03831569
3: VEGF_A_UP.V1_DN 167 Down 0.04466353
rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Group == 3])
treatment <- rnaseq.subset$meta[,factor(Treatment.arm,levels=c("vehicle","GDC0032"))]
d <- model.matrix(~ 0 + treatment)
d
treatmentvehicle treatmentGDC0032
1 1 0
2 1 0
3 1 0
4 1 0
5 1 0
6 0 1
7 0 1
8 0 1
9 0 1
10 0 1
11 0 1
12 0 1
13 0 1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment
[1] "contr.treatment"
cont.matrix <- matrix(c(-1,1), dim(d)[2], 1)
colnames(cont.matrix) <- c("GDC.minus.control")
rownames(cont.matrix) <- colnames(d)
cont.matrix
GDC.minus.control
treatmentvehicle -1
treatmentGDC0032 1
#Perform DE
DE.data <- DE.run(rnaseq.subset$human,d,5,cont.matrix = cont.matrix)
#Top DE genes table
DE <- topTable.annotated.DE(DE.data$efit,p.val = 1)
'select()' returned 1:many mapping between keys and columns
if(nrow(DE[adj.P.Val<0.05])!=0){
DE <- DE[adj.P.Val<0.05]
table(DE[,direction])
DE[direction=="up.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
DE[direction=="down.reg",.(direction,ENSEMBL,SYMBOL,adj.P.Val,logFC)]
#Pheatmap w/ batch correction
i <- match(DE$ENSEMBL,rownames(DE.data$v$E))
y <- removeBatchEffect(DE.data$v$E[i,],design = model.matrix(~drug), batch= pdx)
labels.col <- rep("",nrow(rnaseq$meta))
labels.row <- DE$SYMBOL
ann.col <- as.data.frame(rnaseq$meta[,.(pdx,drug)])
rownames(ann.col) <- rnaseq$meta[,Sample.name]
pheatmap(y,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
} else {
print("No DE genes found")
head(DE,10)
}
[1] "No DE genes found"
ENSEMBL logFC AveExpr t P.Value adj.P.Val B direction ENTREZID SYMBOL GENENAME
1: ENSG00000163995 0.6837759 3.410971 5.496805 3.477881e-05 0.5074877 0.17215334 up.reg 84448 ABLIM2 actin binding LIM protein family member 2
2: ENSG00000111300 0.4156489 5.496378 4.531102 2.726366e-04 0.5074877 0.10500080 up.reg 80018 NAA25 N(alpha)-acetyltransferase 25, NatB auxiliary subunit
3: ENSG00000102753 0.3786802 6.106180 4.424698 3.438274e-04 0.5074877 0.03975911 up.reg 3839 KPNA3 karyopherin subunit alpha 3
4: ENSG00000253729 0.4086519 6.749566 4.330204 4.227513e-04 0.5074877 -0.05735691 up.reg 5591 PRKDC protein kinase, DNA-activated, catalytic polypeptide
5: ENSG00000204389 -0.5307251 7.717875 -4.208828 5.516668e-04 0.5074877 -0.22188442 down.reg 3303 HSPA1A heat shock protein family A (Hsp70) member 1A
6: ENSG00000139055 -0.4244053 5.203128 -4.308408 4.434221e-04 0.5074877 -0.26761439 down.reg 121506 ERP27 endoplasmic reticulum protein 27
7: ENSG00000204922 -0.4767766 4.225305 -4.500622 2.913477e-04 0.5074877 -0.27855530 down.reg 790955 UQCC3 ubiquinol-cytochrome c reductase complex assembly factor 3
8: ENSG00000072501 0.2871578 7.901466 4.138382 6.440253e-04 0.5074877 -0.34573791 up.reg 8243 SMC1A structural maintenance of chromosomes 1A
9: ENSG00000163346 -0.3361114 6.606192 -4.118766 6.724169e-04 0.5074877 -0.41793193 down.reg 57326 PBXIP1 PBX homeobox interacting protein 1
10: ENSG00000234444 0.4092223 4.950297 4.180120 5.875733e-04 0.5074877 -0.57590263 up.reg 728927 ZNF736 zinc finger protein 736
#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- merge(y,rnaseq.subset$meta)
ggplot(y) + aes(x=SYMBOL,y=Expression,color=Treatment.arm) + geom_boxplot() + stat_compare_means(aes(label=paste0("p.adj = ",..p.adj..)),method="t.test") + theme_tufte() + theme(legend.position = "top")
## Heatmap
labels.col <- rnaseq.subset$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq.subset$meta[,.(Treatment.arm)])
rownames(ann.col) <- rnaseq.subset$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
#GSEA - c6
gsea <- GSEA.run(DE.data$v,d,cont.matrix[,1],Hs.c6,0.05)
'select()' returned 1:many mapping between keys and columns
gsea[Direction=="Up",.(rn,NGenes,Direction,FDR)]
rn NGenes Direction FDR
1: VEGF_A_UP.V1_DN 167 Up 2.498761e-05
2: ERB2_UP.V1_DN 186 Up 4.565218e-05
3: TBK1.DF_DN 267 Up 6.914354e-05
4: PIGF_UP.V1_UP 156 Up 5.032092e-04
5: GCNP_SHH_UP_LATE.V1_UP 157 Up 1.941501e-03
6: E2F1_UP.V1_UP 171 Up 4.030628e-03
7: GCNP_SHH_UP_EARLY.V1_UP 146 Up 4.030628e-03
8: RPS14_DN.V1_DN 154 Up 4.522622e-03
9: LTE2_UP.V1_DN 177 Up 6.406573e-03
10: CSR_EARLY_UP.V1_UP 145 Up 7.089887e-03
11: BCAT_BILD_ET_AL_DN 38 Up 1.263871e-02
12: PRC2_EZH2_UP.V1_UP 132 Up 1.722478e-02
gsea[Direction=="Down",.(rn,NGenes,Direction,FDR)]
rn NGenes Direction FDR
1: ERB2_UP.V1_UP 174 Down 0.007089887
2: EIF4E_UP 78 Down 0.027986355
rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Treatment.arm=="vehicle"])
treatment.group <- rnaseq.subset$meta[,Group]
d <- model.matrix(~ 0 + treatment.group)
d
treatment.group1 treatment.group2 treatment.group3
1 1 0 0
2 1 0 0
3 1 0 0
4 1 0 0
5 0 0 1
6 0 0 1
7 0 0 1
8 0 0 1
9 0 0 1
10 0 1 0
11 0 1 0
12 0 1 0
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment.group
[1] "contr.treatment"
cont.matrix <- matrix(c(-1,1,0,-1,0,1,0,-1,1), dim(d)[2], 3)
colnames(cont.matrix) <- c("Long.minus.short","Resistant.minus.short","Resistant.minus.long")
rownames(cont.matrix) <- colnames(d)
cont.matrix
Long.minus.short Resistant.minus.short Resistant.minus.long
treatment.group1 -1 -1 0
treatment.group2 1 0 -1
treatment.group3 0 1 1
#Perform DE
DE.data <- DE.run(rnaseq.subset$human,d,3,cont.matrix = cont.matrix)
dt <- decideTests(DE.data$efit)
summary(dt)
Long.minus.short Resistant.minus.short Resistant.minus.long
Down 2243 2816 0
NotSig 12186 11123 16786
Up 2357 2847 0
dt.common <- dt[rowSums(abs(dt))!=0,]
vennDiagram(dt.common)
#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- merge(y,rnaseq.subset$meta)
ggplot(y) + aes(x=SYMBOL,y=Expression,color=Group) + geom_boxplot() + theme_tufte() + theme(legend.position = "top")
## Heatmap
labels.col <- rnaseq.subset$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq.subset$meta[,.(Group)])
rownames(ann.col) <- rnaseq.subset$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
rnaseq.subset <- subset.PDX(rnaseq, rnaseq$meta[,Treatment.arm=="GDC0032"])
treatment.group <- rnaseq.subset$meta[,Group]
d <- model.matrix(~ 0 + treatment.group)
d
treatment.group1 treatment.group2 treatment.group3
1 1 0 0
2 1 0 0
3 1 0 0
4 1 0 0
5 0 0 1
6 0 0 1
7 0 0 1
8 0 0 1
9 0 0 1
10 0 0 1
11 0 0 1
12 0 0 1
13 0 1 0
14 0 1 0
15 0 1 0
16 0 1 0
attr(,"assign")
[1] 1 1 1
attr(,"contrasts")
attr(,"contrasts")$treatment.group
[1] "contr.treatment"
cont.matrix <- matrix(c(-1,1,0,-1,0,1,0,-1,1), dim(d)[2], 3)
colnames(cont.matrix) <- c("Long.minus.short","Resistant.minus.short","Resistant.minus.long")
rownames(cont.matrix) <- colnames(d)
cont.matrix
Long.minus.short Resistant.minus.short Resistant.minus.long
treatment.group1 -1 -1 0
treatment.group2 1 0 -1
treatment.group3 0 1 1
#Perform DE
DE.data <- DE.run(rnaseq.subset$human,d,4,cont.matrix = cont.matrix)
dt <- decideTests(DE.data$efit)
summary(dt)
Long.minus.short Resistant.minus.short Resistant.minus.long
Down 78 1791 514
NotSig 16695 13186 16005
Up 15 1811 269
dt.common <- dt[rowSums(abs(dt))!=0,]
vennDiagram(dt.common)
#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- merge(y,rnaseq.subset$meta)
ggplot(y) + aes(x=SYMBOL,y=Expression,color=Group) + geom_boxplot() + theme_tufte() + theme(legend.position = "top")
## Heatmap
labels.col <- rnaseq.subset$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq.subset$meta[,.(Group)])
rownames(ann.col) <- rnaseq.subset$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="row")
treatment.group <- rnaseq$meta[,Group]
treatment <- rnaseq$meta[,factor(Treatment.arm,levels=c("vehicle","GDC0032"))]
d <- model.matrix(~ 0 + treatment.group * treatment)
d
treatment.group1 treatment.group2 treatment.group3 treatmentGDC0032 treatment.group2:treatmentGDC0032 treatment.group3:treatmentGDC0032
1 1 0 0 0 0 0
2 1 0 0 0 0 0
3 1 0 0 0 0 0
4 1 0 0 0 0 0
5 1 0 0 1 0 0
6 1 0 0 1 0 0
7 1 0 0 1 0 0
8 1 0 0 1 0 0
9 0 0 1 0 0 0
10 0 0 1 0 0 0
11 0 0 1 0 0 0
12 0 0 1 0 0 0
13 0 0 1 0 0 0
14 0 0 1 1 0 1
15 0 0 1 1 0 1
16 0 0 1 1 0 1
17 0 0 1 1 0 1
18 0 0 1 1 0 1
19 0 0 1 1 0 1
20 0 0 1 1 0 1
21 0 0 1 1 0 1
22 0 1 0 0 0 0
23 0 1 0 0 0 0
24 0 1 0 0 0 0
25 0 1 0 1 1 0
26 0 1 0 1 1 0
27 0 1 0 1 1 0
28 0 1 0 1 1 0
attr(,"assign")
[1] 1 1 1 2 3 3
attr(,"contrasts")
attr(,"contrasts")$treatment.group
[1] "contr.treatment"
attr(,"contrasts")$treatment
[1] "contr.treatment"
DE.data <- DE.run(rnaseq$human,d,coefs = 1:6)
dt <- decideTests(DE.data$efit)
summary(dt)
treatment.group1 treatment.group2 treatment.group3 treatmentGDC0032 treatment.group2:treatmentGDC0032 treatment.group3:treatmentGDC0032
Down 450 185 519 0 1 1
NotSig 2351 2773 2253 17381 17380 17380
Up 14580 14423 14609 0 0 0
#Genes of Interest
x <- DE.data$v$E[match(GoI$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- GoI$SYMBOL
## Boxplot
y <- melt(x)
names(y) <- c("SYMBOL","fname","Expression")
y <- data.table(merge(y,rnaseq$meta))
y[,Treatment.arm:=factor(Treatment.arm,levels=c("vehicle","GDC0032"))]
levels(y$Treatment.arm) <- c("Vehicle","GDC0032")
ggplot(y) + aes(x=Treatment.arm,y=Expression,color=Group) + geom_boxplot() + theme_tufte() + theme(legend.position = "top") + facet_wrap(~SYMBOL)
ggplot(y) + aes(x=Group,y=Expression,color=Treatment.arm) + geom_boxplot() + theme_tufte() + theme(legend.position = "top") + facet_wrap(~SYMBOL)
## Heatmap
labels.col <- rnaseq$meta[,Barcode]
labels.row <- GoI$SYMBOL
ann.col <- as.data.frame(rnaseq$meta[,.(Group,Treatment.arm)])
rownames(ann.col) <- rnaseq$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="none")
# Pathways of Interest
plts <- lapply(PoI, function(p){
x <- DE.data$v$E[match(p$ENSEMBL,rownames(DE.data$v$E)),]
rownames(x) <- p$SYMBOL
x <- x[rowSums(is.na(x)) == 0,]
labels.col <- rnaseq$meta[,Barcode]
labels.row <- rownames(x)
ann.col <- as.data.frame(rnaseq$meta[,.(Group,Treatment.arm)])
rownames(ann.col) <- rnaseq$meta[,fname]
pheatmap(x,labels_col=labels.col,labels_row=labels.row,annotation_col = ann.col,scale="none")
})
ggplot(y[y$SYMBOL=="LDHA",]) + aes(x=Group,y=Expression,fill=Treatment.arm) + geom_boxplot() + stat_compare_means(label="p.format",method="t.test",method.args = list(alternative = "less"),label.y=9.7) + theme_tufte(12) + theme(legend.position = "top",legend.title=element_blank()) + ylab("LDHA Expression")
ggplot(y[y$SYMBOL=="FOXM1",]) + aes(x=Group,y=Expression,fill=Treatment.arm) + geom_boxplot() + stat_compare_means(label="p.format",method="t.test",method.args = list(alternative = "less"),label.y=9.7) + theme_tufte(12) + theme(legend.position = "top",legend.title=element_blank()) + ylab("FOXM1 Expression")
my_comparisons <- list(c("1","2"),c("2","3"),c("1","3"))
ggplot(y[y$SYMBOL=="FOXM1" & Treatment.arm=="GDC0032",]) + aes(x=Group,y=Expression) + geom_boxplot() + stat_compare_means(label.y=8.5) + stat_compare_means(comparisons = my_comparisons,method.args = list(alternative = "less")) + theme_tufte(12) + theme(legend.position = "top") + ylab("FOXM1 Expression") + xlab("GDC-0032 Treated Group")
z <- merge(data.table(t(x),keep.rownames = "fname"),rnaseq$meta)
ggpairs(z,columns = 2:6)